Calibración de probabilidades

Estimación de default mediante modelos de machine learning

Karina Bartolomé

Organizadora: Natalia Salaberry
CIMBAGE (IADCOM) - Facultad Ciencias Económicas (UBA)

2024-05-13

Librerías 
import pandas as pd
import numpy as np
from copy import deepcopy
import seaborn as sns
import matplotlib.pyplot as plt
from collections import Counter
from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.pipeline import Pipeline
from sklearn.compose import ColumnTransformer
from sklearn.compose import make_column_selector
from sklearn.preprocessing import StandardScaler, MinMaxScaler, OneHotEncoder
from sklearn.impute import SimpleImputer
from sklearn.linear_model import LogisticRegression
from sklearn.tree import DecisionTreeClassifier, plot_tree
from sklearn.ensemble import RandomForestClassifier, HistGradientBoostingClassifier
from sklearn.calibration import calibration_curve, CalibratedClassifierCV
from sklearn.metrics import (
    classification_report,
    roc_auc_score,
    roc_curve,
    RocCurveDisplay,
    log_loss,
    recall_score,
    brier_score_loss,
    confusion_matrix,
    ConfusionMatrixDisplay,
)
from IPython.display import display
from custom_functions import (
    plot_distribution, plot_calibration_models, gain_table
)
from great_tables import GT, from_column, style, loc
from matplotlib.colors import ListedColormap, to_rgb

import sklearn
sklearn.set_config(transform_output="pandas")

1 Introducción a machine learning

1.1 ¿Qué es un modelo de clasificación?

Se busca predecir la probabilidad de ocurrencia de un evento a partir de ciertas características observables:

\(P(y=1) = f(X)\)

Siendo y una variable que puede tomar 2 valores: 0 o 1


Ejemplos de clasificación:

▪️Iris: Clasificación de especies de plantas

▪️Titanic: Probabilidad de supervivencia

▪️German Credit: Probabilidad de default

2 Caso: German Credit

2.1 Datos

Código 
PATH_DATA = 'https://raw.githubusercontent.com/ayseceyda/german-credit-gini-analysis/master/gc.csv'
TARGET = "Risk"

df = (pd.read_csv(PATH_DATA)
  .drop('Unnamed: 0', axis=1)
  .assign(Risk = lambda x: np.where(x['Risk']=='good',0,1))
)

Se cuenta con un dataset de 1000 observaciones y 10 variables.

Código 
def map_color(data):
    return (data[TARGET] == 1).map(
        {True: color_verde_claro, False: 'white'}
    )

(GT(df.sample(6, random_state=123)
        .set_index(TARGET)
        .reset_index()
    )
    .fmt_currency(columns="Credit amount")
    .tab_style(
        style=style.fill(color=map_color), locations=loc.body(columns=df.columns.tolist()),
    )
    .tab_style(
        style=style.text(color='red', weight = "bold"), locations=loc.body(TARGET),
    )
    .tab_options(
        column_labels_background_color=color_verde,
        table_font_names="Times New Roman"
    )
)
Tabla 1: Datos de German Credit (muestra de 6 observaciones)
Risk Age Sex Job Housing Saving accounts Checking account Credit amount Duration Purpose
1 29 male 2 own little little $6,887.00 36 education
1 21 male 2 rent little little $902.00 12 education
0 29 male 1 own moderate $2,333.00 24 furniture/equipment
1 20 female 2 rent little little $2,039.00 18 furniture/equipment
0 35 male 2 own moderate $2,728.00 15 radio/TV
0 23 male 2 own moderate $1,444.00 15 radio/TV

La variable objetivo (target) es Risk, donde el porcentaje de observaciones de clase 1 (riesgosos) es 30.0%

2.2 Particiones

Partición en dataset de entrenamiento, validación y evaluación.

  • Dataset de entrenamiento –> Ajuste del modelo
  • Dataset de validación –> Calibración del modelo
  • Dataset de evaluación –> Métricas del modelo calibrado
Código 
y = df[TARGET]
X = df.drop([TARGET], axis=1)

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.5, shuffle=True, stratify=y, random_state=42
)

# Lo ideal es utilizar la partición de validación para calibración
X_valid, X_test, y_valid, y_test = train_test_split(
    X_test, y_test, test_size=0.5, shuffle=True, stratify=y_test, random_state=42
)

print(f"N observaciones en entrenamiento: {X_train.shape[0]}")
print(f"N observaciones en validación: {X_valid.shape[0]}")
print(f"N observaciones en evaluación: {X_test.shape[0]}")
N observaciones en entrenamiento: 500
N observaciones en validación: 250
N observaciones en evaluación: 250

2.3 Modelo de clasificación simple

Se consideran 2 variables para ajustar un modelo de arbol de decisión con máxima profundidad=2

\[ P(\text{Risk}=1) = f(\text{Credit amount}, \text{Age}) \]


Código 
subset_cols = ['Age', 'Credit amount']

(GT(pd.concat([y_train, X_train], axis=1)
        .sample(6, random_state=123)
        .set_index(TARGET)
        [subset_cols]
        .reset_index()
    )
    .fmt_currency(columns="Credit amount")
    .tab_style(
        style=style.fill(color=map_color), locations=loc.body(columns=df.columns.tolist()),
    )
    .tab_style(
        style=style.text(color='red', weight = "bold"), locations=loc.body(TARGET),
    )
    .tab_options(
        column_labels_background_color=color_verde,
        table_font_names="Times New Roman"
    )
)
Tabla 2: Datos de German Credit (2 variables)
Risk Age Credit amount
0 37 $1,287.00
1 27 $8,648.00
1 34 $1,337.00
0 48 $2,134.00
1 48 $1,082.00
0 53 $2,424.00

Ajuste del modelo:

Código 
X_train_subset = X_train[subset_cols].copy()
clf = DecisionTreeClassifier(max_depth=2).fit(X_train_subset, y_train)
clf
DecisionTreeClassifier(max_depth=2)
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.
Código 
colors = [color_verde, 'red']
labels = ['No riesgoso','Riesgoso']
plt.figure(figsize=(6.5,6.5))
arbol = plot_tree(
    clf,
    filled=True,
    impurity=False,
    feature_names=X_train_subset.columns.tolist(),
    class_names=labels,
    fontsize=8
)
for i, impurity, value in zip(arbol, clf.tree_.impurity, clf.tree_.value):
    # let the max value decide the color; whiten the color depending on impurity (gini)
    r, g, b = to_rgb(colors[np.argmax(value)])
    f = impurity * 2 # for N colors: f = impurity * N/(N-1) if N>1 else 0
    i.get_bbox_patch().set_facecolor((f + (1-f)*r, f + (1-f)*g, f + (1-f)*b))
    i.get_bbox_patch().set_edgecolor('black')
plt.show()
Figura 1: Modelo de árbol de decisión

2.4 Evaluación de un modelo de clasificación

Se utiliza el método clf.predict()para predecir la clase (0,1) según el modelo (clf). A partir de la clase predicha se obtiene la matriz de confusión:

Ver código 
preds = pd.DataFrame({'Valor observado':y_test,'Valor predicho':clf.predict(X_test[subset_cols])})
display(GT(preds
    .groupby(['Valor observado','Valor predicho'], as_index=False).size()
    .pivot(index='Valor observado', columns='Valor predicho', values='size')
    .rename({0:'0',1:'1'}, axis=1)
    .reset_index()
    )
    .tab_spanner('Valor predicho', columns = ['0','1'])
    .data_color(
        columns=['0','1'],
        domain=[0,X.shape[0]],
        palette=[color_verde_claro, 'red'],
        na_color="white",
    )
)
Tabla 3: Matriz de confusión
Valor observado Valor predicho
0 1
0 172 3
1 71 4

Además de predecir una clase, sklearn permite predecir una “probabilidad” con clf.predict_proba():

Ver código 
y_pred_proba = clf.predict_proba(X_test[subset_cols])[:, 1]
preds = pd.DataFrame({
    'y_obs':y_test,
    'predict_proba':y_pred_proba
})
GT(preds.sample(2, random_state=42)).fmt_number('predict_proba')
Tabla 4: clf.predict_proba(), muestra de 2 observaciones
y_obs predict_proba
0 0.25
1 0.33 

ROC AUC = 0.64, Log loss = 0.59, Brier loss = 0.2
Ver código 
plot_distribution(data=preds, pred_column='predict_proba', class_0_color=color_verde)
Figura 2: Distribución de “probabilidad” predicha

2.5 Preprocesamiento

Se construye un pipeline de preprocesamiento de variables, diferenciando el procesamiento de variables categóricas y numéricas

Ver código 
numeric_transformer = Pipeline([
    ('impute', SimpleImputer(strategy='median')),
    ('scaler', MinMaxScaler())
])

categorical_transformer = Pipeline([
    ('ohe',OneHotEncoder(
        min_frequency=0.05,
        handle_unknown='infrequent_if_exist',
        sparse_output=False)
    )
])

preproc = ColumnTransformer([
    ('num', numeric_transformer,
      make_column_selector(dtype_include=['float','int'])),
    ('cat', categorical_transformer,
      make_column_selector(dtype_include=['object','category']))
], verbose_feature_names_out=False)
ColumnTransformer(transformers=[('num',
                                 Pipeline(steps=[('impute',
                                                  SimpleImputer(strategy='median')),
                                                 ('scaler', MinMaxScaler())]),
                                 <sklearn.compose._column_transformer.make_column_selector object at 0x14a96f6a0>),
                                ('cat',
                                 Pipeline(steps=[('ohe',
                                                  OneHotEncoder(handle_unknown='infrequent_if_exist',
                                                                min_frequency=0.05,
                                                                sparse_output=False))]),
                                 <sklearn.compose._column_transformer.make_column_selector object at 0x14a96faf0>)],
                  verbose_feature_names_out=False)
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

2.6 Modelado

Pipeline(steps=[('preproc',
                 ColumnTransformer(transformers=[('num',
                                                  Pipeline(steps=[('impute',
                                                                   SimpleImputer(strategy='median')),
                                                                  ('scaler',
                                                                   MinMaxScaler())]),
                                                  <sklearn.compose._column_transformer.make_column_selector object at 0x14a96f6a0>),
                                                 ('cat',
                                                  Pipeline(steps=[('ohe',
                                                                   OneHotEncoder(handle_unknown='infrequent_if_exist',
                                                                                 min_frequency=0.05,
                                                                                 sparse_output=False))]),
                                                  <sklearn.compose._column_transformer.make_column_selector object at 0x14a96faf0>)],
                                   verbose_feature_names_out=False)),
                ('model',
                 HistGradientBoostingClassifier(class_weight={0: 1, 1: 1},
                                                max_depth=4, max_iter=1000,
                                                random_state=42))])
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

2.7 Evaluación

              precision    recall  f1-score   support

           0       0.80      0.79      0.80       175
           1       0.53      0.53      0.53        75

    accuracy                           0.72       250
   macro avg       0.66      0.66      0.66       250
weighted avg       0.72      0.72      0.72       250


ROC AUC = 0.74
Log loss = 0.93
Brier loss = 0.22

2.8 Problema

Planteo de problema organizacional

3 Calibración de probabilidades

3.1 Calibración de probabilidades

Objetivo:

Se busca encontrar una función que ajuste la relación entre los scores predichos por el modelo (predict_proba) y las probabilidades reales

\(P(y_i=1)= f(z)\)

Siendo: \(P(y_{i}=1)\) la probabilidad calibrada para el individuo \(i\) y \(z_{i}\) el output del modelo no calibrado (score)


Un clasificador binario bien calibrado debería clasificar de forma tal que para las observaciones que predice un score (predict_proba) = 0.8, se espera que un 80% de las observaciones correspondan a la clase positiva (1).


Referencias:


Métodos de calibración

Calibración Sigmoide (Platt scaling)

Calibración Isotónica

Calibración Beta (Nuevo)

Calibración Spline (Nuevo)

3.2 Calibración sigmoide (Platt scaling)

Este método asume una relación logística entre los scores (z) y la probabilidad real (p):

\(log(\frac{p}{1-p})=\alpha+\beta(z_{i})\)

\(P(y_{i}=1) = \frac{1}{1+(e^{-(α+β(z_{i}))})}\)

Siendo: \(y_{i}\) el valor observado para el individuo \(i\) y \(z_{i}\) el output del modelo no calibrado


Se estiman 2 parámetros (\(\alpha\) y \(\beta\)), como en una regresión logística.

  • Requiere pocos datos

  • Es útil cuando el modelo no calibrado tiene errores similares en predicción de valores bajos y altos


Referencias:

  • Platt, John. (2000). Probabilistic Outputs for Support Vector Machines and Comparisons to Regularized Likelihood Methods. Adv. Large Margin Classif, Volume 10 (pp. 61-74).

3.3 Calibración sigmoide (Platt scaling) (Cont.)

Ver código 
cal_data = plot_calibration_models(
    y_obs=preds["y_obs"],
    y_pred_prob=preds["pred_prob"],
    y_pred_prob_cal=preds["pred_sigmoid"],
    bins=7,
    strategy="uniform",
    figsize=(4,4),
    display_bins=True
)
plt.legend(loc='center left', bbox_to_anchor=(0, -0.3))
plt.show()

# display(GT(cal_data).fmt_number(['prob_true','prob_pred']))

plt.figure(figsize=(4,4))
sns.scatterplot(
    x="pred_prob", y="pred_sigmoid", data=preds,
    color="blue", alpha=0.7, label="Calibración sklearn",
)
sns.scatterplot(
    x="pred_prob", y="y_obs", data=preds,
    alpha=0.3, color="black", label="Valores observados",
)
plt.axline([0, 0], [1, 1], label="Calibración perfecta", c="red")
plt.ylabel("Probabilidad calibrada")
plt.xlabel("Probabilidad no calibrada")
plt.legend(loc='center left', bbox_to_anchor=(0, -0.3))
plt.ylim([-0.01, 1.01]);
Figura 3: Calibración (bins)
Figura 4: Calibración (ajuste)

3.4 Calibración isotónica

CalibratedClassifierCV(cv='prefit',
                       estimator=Pipeline(steps=[('preproc',
                                                  ColumnTransformer(transformers=[('num',
                                                                                   Pipeline(steps=[('impute',
                                                                                                    SimpleImputer(strategy='median')),
                                                                                                   ('scaler',
                                                                                                    MinMaxScaler())]),
                                                                                   <sklearn.compose._column_transformer.make_column_selector object at 0x14a96f6a0>),
                                                                                  ('cat',
                                                                                   Pipeline(steps=[('ohe',
                                                                                                    OneHotEncoder(handle_unknown='infrequent_if_exist',
                                                                                                                  min_frequency=0.05,
                                                                                                                  sparse_output=False))]),
                                                                                   <sklearn.compose._column_transformer.make_column_selector object at 0x14a96faf0>)],
                                                                    verbose_feature_names_out=False)),
                                                 ('model',
                                                  HistGradientBoostingClassifier(class_weight={0: 1,
                                                                                               1: 1},
                                                                                 max_depth=4,
                                                                                 max_iter=1000,
                                                                                 random_state=42))]),
                       method='isotonic')
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.
% observado: 0.3
Prob promedio: 0.308
Prob promedio (calibración sigmoide): 0.314
Prob promedio (calibración isotónica): 0.329

3.5 Comparación de modelos

ROC AUC modelo=0.744
ROC AUC modelo calibrado=0.744
Log loss modelo=0.929
Log loss modelo calibrado=0.531
Brier modelo=0.223
Brier modelo calibrado=0.177
Recall modelo=0.533
Recall modelo calibrado=0.400

Se visualizan las distribuciones de probabilidad predicha según el modelo calibrado y no calibrado.

Ver código 
plot_distribution(
    data=preds, pred_column='pred_prob'
)
plot_distribution(
    data=preds, pred_column='pred_sigmoid'
)

Modelo no calibrado

Modelo calibrado (sigmoide)

3.6 Visualización

3.7 Tablas de ganancias

Se busca un modelo en donde la probabilidad promedio de cada bin se corresponda con el % de clase positiva observado en ese bin según las predicciones del modelo.

Esto es útil para la toma de decisiones, ya que permite establecer punto de cortes diferenciales en función de la aversión en riesgo de la entidad.

Por ejemplo:

  • En un modelo de scoring crediticio, otorgarle créditos a N individuos con probabilidad en cierto intervalo tiene un riesgo asociado (mora esperada)
  • En un modelo de churn (abandono) de clientes, otorgarle una promoción a individuos de cierto intervalo de probabilidad de abandono tiene un costo asociado (descuento para individuos que no abandonarían)
  • Entre otros.

3.8 Tabla de ganancias (Cont.)

Código 
gain_table(preds=preds, bins=5)
Tabla 5: Tabla de ganancias (modelo no calibrado)
Bin N N_clase1 avg_prob perc_clase1 diff N_acum N_clase1_acum TPR
(-0.10%, 0.12%] 50 4 0.03% 8.00% 0.080 50 4 8.00%
(0.12%, 1.75%] 50 7 0.65% 14.00% 0.134 100 11 11.00%
(1.75%, 22.80%] 50 14 7.58% 28.00% 0.204 150 25 16.67%
(22.80%, 83.70%] 50 20 50.60% 40.00% 0.106 200 45 22.50%
(83.70%, 100.00%] 50 30 95.09% 60.00% 0.351 250 75 30.00%
Código 
gain_table(preds=preds, pred_column='pred_sigmoid', bins=5)
Tabla 6: Tabla de ganancias (calibración sigmoide)
Bin N N_clase1 avg_prob perc_clase1 diff N_acum N_clase1_acum TPR
(1.20%, 14.10%] 50 4 8.69% 8.00% 0.007 50 4 8.00%
(14.10%, 22.80%] 50 7 18.47% 14.00% 0.045 100 11 11.00%
(22.80%, 35.00%] 50 14 28.12% 28.00% 0.001 150 25 16.67%
(35.00%, 49.90%] 50 20 41.46% 40.00% 0.015 200 45 22.50%
(49.90%, 80.70%] 50 30 60.25% 60.00% 0.003 250 75 30.00%

Facultad de Ciencias Económicas - UBA |